Kinetic Monte Carlo Simulations of a Model for Heat-assisted Magnetization Reversal 

in Ultrathin Films 
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To develop practically useful systems for ultra-high-density information recording with densi- 
ties above terabits/cm 2 , it is necessary to simultaneously achieve high thermal stability at room 
temperature and high recording rates. One method that has been proposed to reach this goal is 
heat-assisted magnetization reversal (HAMR). In this method, the magnetic orientation is assigned 
to a high-coercivity material by temporarily reducing the coercivity during the writing process 
through localized heating. Here we present kinetic Monte Carlo simulations of a model of HAMR 
for ultrathin films, in which the temperature in the central part of the film is momentarily in- 
creased above the critical temperature, for example by a laser pulse. We observe that the speed-up 
achieved by this method, relative to the switching time at a constant, subcritical temperature, is 
optimal for an intermediate strength of the writing field. This effect is explained using the theory 
of nucleation-induced magnetization switching in finite systems. Our results should be particularly 
relevant to recording media with strong perpendicular anisotropy, such as ultrathin Co/Pt or Co/Pd 
multilayers. 

PACS numbers: 75.60.Jk 75.70.Cn 05.70.Ln 64.60.Dc 



I. INTRODUCTION 

One of the most important factors supporting progress 
in the miniaturization of computers and other elec- 
tronic devices is the continued exponential increase in 
the density of data storage.— Currently, designs are 
being considered for magnetic recording devices that 
have areal data densities of the order of terabits/cm 2 
- several orders of magnitude more than only a decade 
ago. At such densities, the size of the recording bit 
approaches the superparamagnetic limit, where ther- 
mal fluctuations seriously degrade the stability of the 
magnetizationi 2 * 3 - However, current industry standards 
demand that bits should retain 95% of their magneti- 
zation over a period of ten yearsi Furthermore, sub- 
nanosecond magnetization-switching times are required 
to achieve acceptable read/write rates. 

One suggested method to fulfill these requirements is 
to use ultrathin, perpendicularly magnetized films of very 
high-coercivity materials, such as FePt (coercive field 
about 50 kOe), or single-particle bits that are expected 
to have even higher coercivities i However, such high co- 
ercive fields at room temperature are beyond what is 
achievable by modern write heads, which are limited to 
about 17 kOe.— A method suggested to overcome this 
problem is to exploit the temperature dependence of the 
coercivity through heat-assisted magnetization reversal 



or HAMR (aka. thermally assisted magnetization rever- 
sal or TAMR)jiiir— This is accomplished by increasing 
the temperature of the recording area to a value close 
to, or above, the Curie temperature of the medium via 
a localized heat source, such as a laser Due to 

the temperature dependence of the coercivity, the mag- 
nitude of the required switching field is lowered at the 
elevated temperature, relaxing the requirements for the 
write head. An important consideration for the imple- 
mentation of the HAMR technique is to keep the heat 
input as low and as tightly focused as possible, limiting 
energy transfer to neighboring recording bits. In order to 
reach the desired high data densities, the laser spot must 
have a diameter less than 50 nm, much smaller than the 
wavelength. This can be achieved using near-field optics, 
a technology which currently is the objective of vigorous 
research and development &12r— 

Despite their simplicity, two-dimensional kinetic Ising 
models have been shown to be useful for studying 
magnetization switching in ultrathin films with strong 
anisotropy. 3 Theoretical^ 3 - and experimental^ work has 
shown that the equilibrium phase transition in such films 
belongs to the universality class of the two-dimensional 
Ising model. The dynamics of magnetization switching 
in ultrathin, perpendicularly magnetized films has been 
studied using magneto-optical microscopies in combina- 
tion with Monte Carlo simulations of Ising-like models 



by, among others, Kirilyuk et alJ^ and Robb et alJ£ Sys- 
tems that have been found to have strong Ising charac- 
ter include Fe sesquilayersi 4 - and ultrathin films of Co,— 
Co/Pd^i and Co/Pt i 16 ' 18 The strong anisotropy in such 
systems limits the effects of transverse spin dynamics and 
ensures that local spin reversals are thermally activated. 
The extreme thinness of the films strongly reduce the 
demagnetization effects to which films with out-of plane 
magnetization are otherwise subject] 13 ' 14 ' 16 For detailed 
reviews of experimental and simulational studies of mag- 
netization switching in ultrathin films with perpendicular 
magnetization, see Refs. [l8lfl9l . 

In the present paper we use a two-dimensional Ising fcr- 
romagnet to model the HAMR process by kinetic Monte 
Carlo (MC) simulation, demonstrating enhanced nucle- 
ation of the switched magnetization state in the heated 
area. For simplicity and computational economy, we en- 
visage an experimental setup slightly different from oth- 
ers previously reported in the literature It most 
closely resembles the optical-dominant setup shown in 
Fig. 1(b) of Ref. 1. The recording medium is placed 
in a constant write field that is too weak to cause sig- 
nificant switching on an acceptable time scale, and it 
is heated at its center by a transient heat pulse. At a 
fixed superheating temperature we show that the rela- 
tive speed-up of the magnetization switching, compared 
to the constant-temperature case, depends nonmonoton- 
ically on the magnitude of the applied field. This relative 
speed-up shows a pronounced maximum at an interme- 
diate value of the applied field. We give a physical ex- 
planation for this effect, based on the nucleation theory 
of magnetization switching in finite-sized systems! 3 ' 20 ' 21 
As magnetization switching is a special case of the de- 
cay of a metastable phase (i.e., the medium in its state 
of magnetization opposite to the applied field) ) 21 ' 22 this 
analysis is of general physical interest beyond the specific 
technological application discussed here. 

The rest of this paper is organized as follows. Our 
model and methods are described in Sec. [Til the numer- 
ical results are described and explained in Sec. IIII1 and 
our conclusions are stated in Sec. HV] 



II. MODEL AND METHODS 

We use a square-lattice, nearest-neighbor Ising ferro- 
magnet with energy given by the Hamiltonian, 

Here, Si = ±1, J > is the strength of the spin inter- 
actions, and the first sum runs over all nearest-neighbor 
pairs. For convenience we hereafter set J = 1. In the 
second term, which represents the Zeeman energy, H is 
proportional to a uniform external magnetic field, and 
the sum runs over all lattice sites. We use a lattice of size 
L 2 = 128 x 128, with periodic boundary conditions. The 



length unit used in this study is the computational lattice 
constant, which should correspond to a few nanometers. 

For simplicity, our model does not include any explicit 
randomness, such as impurities or random interaction 
strengths. As a result, pinning of interfaces for very 
weak applied field, 15,16 as well as heterogeneous nucle- 
ation of spin reversal^ are neglected. We further ex- 
clude demagnetizing effects, which are very weak for ul- 
trathin film s 13 ' 14 ' 16 and thus cause no qualitative changes 
in Monte Carlo simulations of the switching process. 20 

The stochastic spin dynamic is given by the single-spin 
flip Metropolis algorithm with transition probability^ 3 - 

P(s l -> -Si) = min[l, exp(-AE/T)] , (2) 

where AE is the energy change that would result from 
acceptance of the proposed spin flip. The temperature, 
T, is given in energy units (i.e., Boltzmann's constant 
is taken as unity). Updates are attempted for randomly 
chosen spins, and L 2 attempts constitute one MC step 
per spin (MCSS), which is the time unit used in this 
work. (We note that the Metropolis algorithm is not the 
only Monte Carlo dynamics that could be used here. We 
have chosen it because of its simplicity and ubiquity in 
the literature since we do not expect that the inclusion of 
complications such intrinsic barriers to single-spin flips 
would have significant effects at this high temperature 
beyond a renormalization of the overall timescale.) 

Following this algorithm and starting from Sj = — 1 
for all i, we equilibrate the system over 4 x 10 4 MCSS 
at H = and temperature Tq = 0.8T C ss 1.82, where 
T c = 2/ ln(l+\/2) = 2.269... is the exact critical tempera- 
ture for the square-lattice Ising model^ Having achieved 
equilibrium with negative magnetization at zero field, we 
then subject the system to a constant, uniform, posi- 
tive magnetic field, along with a transient heat pulse. 
To simulate the heat pulse, we use a temperature profile 
given by a time-dependent, Gaussian solution of a one- 
dimensional diffusion equation. The profile is centered 
on the mid- line of the Ising lattice, x = 63.5, and each 
spin in the ccth column of the lattice has the temperature 

T ^^ To + - 3T ^ eXP (-4fef))' 

(3) 

Here, 0.3T C is the maximum of the temperature pulse, 
which is attained at t = 0. Therefore, the peak tempera- 
ture is Xo+0.3T c = 1.1T C . The parameter k is the thermal 
diffusivity, which is also set to unity for convenience. The 
time to = a 2 /2k is related to the duration of the heat- 
input process, such that a is the standard deviation that 
governs the width of the temperature profile at t = 0^ 
Here we use a — 6 for all simulations. (Equation [3] most 
likely underestimates the speed of decay of the temper- 
ature pulse as it ignores heat conduction into the sub- 
strate.) Figure Q] displays the temperature of each col- 
umn at eight times between t = 1 and 500 MCSS. By first 
promoting the center-most lattice sites to temperatures 
above T c before relaxing them back to Tq according to 
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Eq. ([3]), we expect to initiate a magnetization-switching 
event that originates along the center line of the lattice 
and propagates outward. After the completion of this 
switching process, almost all spins will be oriented up, 
S{ = +1. We dchne the switching time t s as the time 
until the system first reaches a magnetization per spin, 

m = j^J2 s *> ( 4 ) 

i 

of zero or greater. 

III. RESULTS 

We first performed a preliminary study to confirm that 
magnetization switching can be induced by the temper- 
ature profile, given the parameters used in Eq. ([3]). For 
this purpose, we inspected snapshots of the system dur- 
ing a single run at H = 0.2. In Fig. [5] we display the 
configuration of the system at six times between t = 1 
and 125 MCSS during this run. As expected, the switch- 
ing begins near the center line of the system, where the 
temperature is above critical, and propagates outward. 
We note a strong similarity of the simulated magnetiza- 
tion configurations to experimental images of ultrathin, 
strongly anisotropic films undergoing magnetization re- 
versal, such as Figs. 3, 4, and 8 of Ref. [15| and Fig. 2 
of Ref. [l6|. This observation further confirms the ability 
of our simplified model to elucidate generic dynamical 
features of real ultrathin films. 

Having confirmed a switching event at H = 0.2, statis- 
tics were accumulated for 200 simulations at H = 0.2 and 
also at fifteen weaker fields down to H = 0.06, as detailed 
in Table|TJ For each field, 100 simulations were performed 
at a constant, uniform temperature of Tq = 0.8T C , and 
100 were performed using the time-dependent tempera- 
ture profile given by Eq. (j3|). For each run, the aver- 
age magnetization for each column at each time step was 
recorded along with the switching time, t s . 

To investigate the effect that the relaxing temperature 
profile has on each column of the Ising lattice, we plotted 
the average magnetization per spin against the column 
number. In Fig. [3] we show this average magnetization for 
H = 0.2, 0.08, and 0.06. The plots on the left [Fig.^a), 

(c) , and (e)] result from the 100 runs with the relaxing 
temperature profile, and the ones on the right [Fig.[3Jb), 

(d) , and (f)] from the 100 runs at the constant, uniform 
temperature of To. The plots at H = 0.2 [(a) and (b)] 
show the average magnetization per spin at eight different 
times between t = 1 and 300 MCSS. The plots at H = 
0.08 [(c) and (d)] show the average magnetization per 
spin at ten different times between t = 1 and 5500 MCSS. 
Finally, the plots at H = 0.06 [(e) and (f)] show the 
average magnetization per spin at nine different times 
between t = 1 and 25000 MCSS. (For a full listing of the 
times, see the figure caption.) 

Again comparing the results with a relaxing tempera- 
ture profile to those realized at constant, uniform tem- 



perature, in Fig. [4] we show cumulative probability distri- 
butions for the switching times for fields H — 0.2, 0.15, 
0.08, 0.0725, 0.065, and 0.06. The black "stairs" are the 
cumulative distributions for the switching times in the 
100 runs with the relaxing temperature profile (hereafter 
referred to as t s ). The gray (red online) stairs are the cu- 
mulative distributions for the switching times in the 100 
runs at constant, uniform temperature (hereafter referred 
to as t c ). 

Table 1 lists the median switching times for both the 
100 runs with the relaxing temperature profile (t s ) and 
the 100 runs at constant, uniform temperature To (^c) for 
each value of H . Also listed are the estimated errors At s 
and At c . The last two columns give the ratio t s /t c and 
the associated error A(t s /t c ). The error At s is defined as 
(i S 2 — i s i)/2, where t S 2 is the switching time with a cu- 
mulative probability of 0.55 and t s i is the switching time 
with a cumulative probability of 0.45, and At c is defined 
analogously. The error in the ratio (t s /t c ) is calculated 
in the standard way as 

The median switching time has the advantage over the 
mean that it can be estimated even when only half of the 
100 simulations switch within the maximum number of 
time steps. This significantly reduces the computational 
requirements, especially for weak fields. 

The ratio (t s /t c ) is plotted vs. H in Fig. [5] The min- 
imum value of this ratio signifies the maximum benefit 
from using the relaxing temperature profile of the HAMR 
method. The corresponding field value, H = 0.0725, is 
the optimal field for this simulation. 

To explain the nonmonotonic shape of the curve rep- 
resenting (t s /t c ) in Fig. [SI it is necessary to understand 
the two most important modes of nucleation-initiated 
magnetization switching in finite-sized systems: mul- 
tidroplet (MD) and single-droplet (SD). (For more de- 
tailed discussions, see Refs. I2lll22h The average time 
between random nucleation events of a growing droplet 
of the equilibrium phase in a d-dimensional system of 
linear size L has the strongly field-dependent form, r n oc 
L d exp[S(T)/(T| J ff d - 1 |)], where S(T) is a measure of the 
free energy associated with the droplet surfaced Once a 
droplet has nucleated, for the weak fields and relatively 
high temperatures studied in this work it grows with a 
near-constant and isotropic radial velocity v g oc \H \/T*& 
As a consequence, the time it would take a newly nu- 
cleated droplet to grow to fill half of a system of vol- 
ume of L d is therefore r g cx L/v g . If r g 3> r n , many 
droplets will nucleate before the first one grows to a size 
comparable to the system, and many droplets will con- 
tribute to the switching process. This is the MD regime, 
which corresponds to moderately strong fields and/or 
large systems^ It is the switching mode shown in Fig. [2] 
for H = 0.2. In the limit of infinitely large systems it is 
identical to the well-known Kolmogorov-Johnson-Mehl- 
Avrami (KJMA) theory of phase transformations^ - — 
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If Tg <C r n , the first droplet to nucleate will switch the 
system magnetization on its own. This is the SD regime, 
which corresponds to weak fields and/or small systems. 21 
It is the switching mode shown in Fig. |B]for H = 0.06. 
The crossover region between the SD and MD regimes is 
known as the Dynamic Spinodal (DSP).— 

One aspect of the MD/SD picture that is particularly 
relevant to the current problem, is the fact that any 
switching event that takes place at a time t < r g cannot 
be accomplished by a single droplet, and thus it must be 
due to the MD mechanism^ For a circular droplet in 
a square L x L system, r g w L/(v / 27rw g ). Using results 
from Ref. [26| (which, like the present model, neglects pin- 
ning effect a 15 ' 16 ), we find that in the range of moderately 
weak fields studied here, at T = 0.8T C v g can be well 
approximated as v g ss 0.75 tanh (H/1.82). The result- 
ing estimates for r g in the simulations (which contain 
no adjustable parameters) are shown as vertical lines in 
Fig. |4jc-f). A kink in the cumulative probability distri- 
bution for the heat-assisted runs is observed at r g , with 
significantly higher slopes in the MD regime on the short- 
time side of r g , than in the SD regime on the long-time 
side. From these figures we see that the optimal field 
value for L = 128, H = 0.0725, corresponds to the situa- 
tion where just above 50% of the heat-assisted switching 
events are caused by the MD mechanism, while essen- 
tially all the constant-temperature switching events are 
SD. This situation is illustrated by the series of snapshots 
in Fig. [7] For significantly larger fields, both protocols 
lead to all MD switching events [Fig. |U[a,b)], while for 
weaker fields, the great majority of the switching events 
are SD for both protocols [Fig.[4je,f)]. In both cases, the 
ratio t s /t c is larger than it is for fields near the optimal 



value [Fig. @Jc,d)]. We have confirmed these conclusions 
by additional simulations for L — 64 and 96 (not shown) . 



IV. CONCLUSIONS 

In this paper we have studied a kinetic Ising model 
of magnetization reversal under the influence of a mo- 
mentary, spatially localized input of energy in the form 
of heat (heat-assisted magnetization reversal, or HAMR). 
Our numerical results indicate that the HAMR technique 
can significantly speed up the magnetization reversal in 
a uniform, applied magnetic field, and we find that this 
speed-up has its optimal value at intermediate values of 
the field. This effect is explained in terms of the MD 
and SD mechanisms of nucleation-initiated magnetiza- 
tion switching in finite systems>21 The two-dimensional 
geometry chosen for this study is particularly appropri- 
ate for thin films. We therefore expect that our predic- 
tions should be experimentally observable for ultrathin 
ferromagnetic films with strong perpendicular anisotropy, 
such as Co/Pd 9 - 17 or Co/Pt 16 - 18 multilayers. 
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(a) H = 0.2, t = 1, m = -0.935 



(b) H = 0.2, t = 10, m = -0.6 



(c) H = 0-2. t = 25, m = -6 791 



(d) H = 0.5. t = 50. m = 569 



TABLE I: Median switching times t s (relaxing temperature 
profile) and t c (constant, uniform temperature) and their ratio 
for the sixteen field values included in Fig. [5] Also given are 
the estimated errors, At s , At c , and A(t a /t c )- 



H 


Median t s 


At s 


Median t c 


At c 


(t s /t c ) A(t s /t c ) 


0.2000 


98.0 


1.0 


126.0 


1.5 


0.778 


0.012 


0.1800 


126.0 


1.5 


166.5 


2.0 


0.757 


0.013 


0.1600 


165.0 


0.5 


225.5 


4.0 


0.732 


0.013 


0.1500 


189.5 


2.5 


263.0 


4.5 


0.721 


0.016 


0.1200 


323.0 


6.0 


482.5 


6.5 


0.669 


0.015 


0.1000 


504.0 


16.0 


881.5 


25.0 


0.572 


0.024 


0.0900 


670.5 


14.5 


1253.0 


92.0 


0.535 


0.041 


0.0800 


1077.0 


43.0 


2015.0 


113.0 


0.534 


0.037 


0.0775 


1148.0 


72.0 


2676.0 


344.5 


0.429 


0.061 


0.0725 


1374.0 


63.5 


4413.0 


354.5 


0.311 


0.029 


0.0700 


2443.0 


274.5 


5470.0 


938.0 


0.447 


0.092 


0.0670 


3232.0 


629.5 


6621.5 


838.5 


0.488 


0.113 


0.0660 


4035.5 


580.5 


7562.5 


733.0 


0.534 


0.093 


0.0650 


6426.5 


1453.0 


9030.0 


1035.5 


0.712 


0.180 


0.0620 


11569.5 


1927.0 


14788.0 


3607.0 


0.782 


0.231 


0.0600 


13808.0 


2479.0 


20851.0 


3435.5 


0.662 


0.161 



(e) H = 02.t = 100. m = -0 023 



(t) H = 0.2, 1 = T2b, m = 0.264 



FIG. 2: Parts (a)-(f) are snapshots of the 128 x 128 Ising 
system at t = 1, 10, 25, 50, 100, and 125 MCSS under in- 
fluence of the time-dependent temperature profile, Eq. 0, 
and a constant, uniform applied field of H — 0.2. Growing 
clusters of the switched phase are first seen to nucleate near 
the center line, where the temperature is highest. However, 
active nucleation is also seen elsewhere in the system. 




FIG. 1: The time dependent Gaussian temperature profile 
used to simulate the decay of a laser heat pulse applied at the 
center line of the Ising lattice. The times plotted are t = 1, 5 
10, 25, 50, 100, 200, 300, and 500 MCSS. The tallest Gaussian 
corresponds to t = 1 MCSS. 
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FIG. 3: The magnetization per spin vs. the column number in 
the lattice, each part [(a)-(f)] averaged over 100 independent 
runs. The plots on the left [(a), (c), and (e)] result from the 
100 runs with the relaxing temperature profile, and the ones 
on the right [(b), (d), and (f)] from the 100 runs at a constant, 
uniform temperature, To = 0.8T C . The plots at H — 0.2 [(a) 
and (b)] show the average magnetization per spin at t = 1, 
10, 25, 50, 100, 150, 200, and 300 MCSS from bottom to 
top. The plots at H — 0.08 [(c) and (d)] show the average 
magnetization per spin at t = 1, 75, 400, 600, 1000, 1300, 
1600, 2000, 3000, and 5500 MCSS from bottom to top. The 
plots at H = 0.06 [(e) and (f)] show the average magnetization 
per spin at t = 1, 500, 1500, 2500, 4000, 7500, 14000, 20000, 
and 25000 MCSS. 
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FIG. 4: Parts (a) - (f) are cumulative probability distribu- 
tions for the switching times with fields H — 0.2, 0.15, 0.08, 
0.0725, 0.065, and 0.06, respectively. The black "stairs" cor- 
respond to the 100 simulations with the relaxing tempera- 
ture profile (switching times, t s ). The gray (red online) stairs 
correspond to the 100 simulations at uniform temperature 
(switching times, t c ). The vertical lines in parts (c) - (f) 
mark the single-droplet growth time r g . Note that the time 
scale increases by more than a factor 100 from (a) to (f). See 
discussion in the text. 
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FIG. 5: The switching-time ratio t s /t c , shown vs. H. The 
minimum value of this ratio signifies the maximum benefit 
from applying the relaxing temperature profile of the HAMR 
method. 
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0.06, t = 500, m = -0.896 




FIG. 6: Parts (a)-(d) are snapshots of the 128 x 128 Ising sys- 
tem at t = 1, 500, 1500, and 2500 MCSS under influence of 
the time-dependent temperature profile, Eq. ([3]), and a con- 
stant, uniform applied field of H = 0.06. In this weak field 
the switching follows the SD mechanism, even in the heated 
region. 



(a) H = 0.0725, t = 1 , m = -0.9466 



(b) H = 0.0725, t = 200, m = -0.865 
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FIG. 7: Parts (a)-(f) are snapshots of the 128 x 128 Ising 
system at t = 1, 200, 400, 600, 800, and 1000 MCSS under 
influence of the time-dependent temperature profile, Eq. ([3|, 
and a constant, uniform applied field of H = 0.075. At this 
intermediate field, multiple growing clusters of the switched 
phase are first seen to nucleate near the center line, where the 
temperature is highest. Without the heat pulse, the switching 
would proceed via the SD mechanism. 
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